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We study the properties of convex functionals which have been proposed for the simulation of charged molec¬ 
ular systems within the Poisson-Boltzmann approximation. We consider the extent to which the functionals 
reproduce the true fluctuations of electrolytes and thus the one-loop correction to mean field theory - includ¬ 
ing the Deby-Hfickel correction to the free energy of ionic solutions. We also compare the functionals for use 
in numerical optimization of a mean field model of a charged polymer and show that different functionals 
have very different stiffnesses leading to substantial differences in accuracy and speed. 


INTRODUCTION 

The sign of electrostatic free energy functionals has 
long interested and puzzled the community of molecu¬ 
lar simulators: A recent paper states The fact that the 
functional cannot be identified with the electrostatic en¬ 
ergy away from the minimum a priori precludes its use 
in a dynamical “on the fly” optimization. . .®. Similar 
statements as to the nature of the free energy in Poisson- 
Boltzmann functionals can be found in other papers 2 -. 

Implicitly most formulations of electrostatic free ener¬ 
gies worlP with the electric field and the potential and 
are associated with the natural potential energy 

Ue = -J~ (!) 

This energy seems to be concave and unbounded below, 
whereas when thermodynamic potentials are written in 
terms of the electric displacement field D we find 

f D 2 

U D = + J — dr (2) 

which is convex. As emphasised in a classic textP the 
two formulations are identical in content, and simply 
linked by a Legendre transform. The importance of the 
correct ensemble in the understanding of electrodynam¬ 
ics and stability of media was particularly emphasised 
in reviews of Kirzhnits and collaborators^ 1 . Inspired 
in particular by this view we proposed a free energy 
functional for th e P oisson-Boltzmann equation which is 
a true minimizei^l, rather than a stationary principle. 
Our hope is that when implemented in practical coded 10 ^ 
these locally formulated convex forms give more stable 
and simpler algorithms. Using a convex free energy it 
is possible to implement an implicit solvent code which 
fully includes the Poisson-Boltzmann free energy in which 
“on the fly” optimization can then be performed using the 
Car-Parrinello method 11 ! In this case the electrodynamic 
field can evolve in parallel to degrees of freedom in the 
moleculeP-^. 

Our work takes the displacement field D as the funda¬ 
mental thermodynamic field, rather than <fi or E. Other 
recent workSHEl takes a different approach to the prob¬ 


lem. By a series of transformations to the original varia¬ 
tional formulations of the Poisson-Boltzmann functional 
(expressed as a function of the electrostatics potential, <jf) 
the authors have found new families of functionals which 
are scalar, convex and local. 

While all these functionals are strictly equivalent at 
their stationary point when simulating at finite temper¬ 
atures differences can occur: The Car-Parrinello method 
thermostats the supplementary variables to zero tem¬ 
perature; in this way we only see the true minimum of 
the energy functional. However in many applications we 
might wish to run simulations with all variables ther- 
mostated to the same, ambient temperature. For electro¬ 
static functionals in dielectric media this samples fluc¬ 
tuations which are equivalent to the so-called thermal 
Casimir interaction 16 ®^!, which is just the zero frequency 
contribution to the free energy in Lifshitz theory. There 
thus arises the question - what happens to the free ener¬ 
gies when a convex Poisson-Boltzmann functional is used 
in a finite temperature simulation? One might hope that 
the functionals also reproduce the correct one-loop cor¬ 
rections to the free energy, which in the case of an elec¬ 
trolyte has an anomalous scaling in p 3 / 2 where p is the 
salt concentration. 

To answer this question we will study the spectra of 
several free energy functionals expanded to quadratic or¬ 
der. The spectra are also crucial in understanding the 
convergence properties and dynamics of simulation algo¬ 
rithms. Accurate (non-implicit) integration in molecular 
dynamics requires a time step which is short compared 
with all dynamic processes under consideration. Algo¬ 
rithms which generate time scales which are very dif¬ 
ferent for different modes are said to be stiff and lead 
to reduced efficiency. Stiffness can also have other, un¬ 
wanted consequences in optimization algorithms - such 
as poor convergence properties. We explore this ques¬ 
tion in the last part of the paper where we study a model 
of free energy minimization in the context of a confined 
polyelectrolyte. 

We start by a reminder of Poisson-Boltzmann theory 
and the corrections to it coming from Gaussian fluctua¬ 
tions at the one-loop level. We then show that the same 
one-loop free energy is found in a complementary dual 
(and convex) formulation of the free energy which can 
be used within a full molecular dynamics or Monte Carlo 
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simulation. We then move on to considerations of coarse 
grained models and test various discretizations of electro¬ 
statics energies for their accuracy and ease of use within 
minimizing principles. 


POISSON-BOLTZMANN THEORY AND ITS ONE-LOOP 
CORRECTION 


In our presentation we will consider the case of a sym¬ 
metric electrolyte, however the identities that we derive 
are independent of the exact microscopic model used. 
The use of a definite physical example leads to consider¬ 
able simplification of notation. 

The mean-field Poisson-Boltzmann functional for a 
symmetric electrolyte is 

= J - 2k B T A cosh f3e(j) + p e <fi^ dr (3) 

where k B T = /3 _1 is the thermal energy, e the ion charge 
<f> the electrostatic potential and A a fugacity. At the 
stationary point of this “free energy” we find 

div egrad <f> — 2Ae sinh ((3ecf>) + p e = 0 (4) 


which is indeed the classic equation for the potential in 
the Poisson-Boltzmann formalism. Linearizing this equa¬ 
tion allows us to define the inverse Debye length from 
k 2 = 2A/3e 2 /e. At low concentrations A can be identified 
with the salt concentration. Since this thermodynamic 
potential is naturally expressed in terms of the thermo¬ 
dynamics fields (j) or E it is concave as noted for eq. ([lj. 
Any attempt to use such a form in Car-Parrinello simu¬ 
lation will lead to numerical instabilities. 

Much recent work has used field theory techniques to 
extend this functional to include fluctuations. In particu¬ 
lar one can sum the fluctuations to the one-loop level 19 20 , 
which gives as a correction to this free energy the func¬ 
tional determinant 


Fioop = —— lo g | - div e(V)grad + 2A/?e 2 cosh {/3e(f>) \ 

( 5 ) 

This expression is often interpreted with a substraction 
scheme - one compares with the free energy of an empty 
box, where A = 0. Even with this subtraction the ex¬ 
pression is formally divergent. However on introducing a 
microscopic cut-o flP the loop free-energy gives both the 
Born self-energy of solvation 


E 


Born 


87rea 


( 6 ) 


as well as the anomalous Debye-Hiickel free energy 


Fdh _ k B Tn 3 
V ~ 12tt 


( 7 ) 


PROPERTIES OF THE LEGENDRE TRANSFORM 


As noted above many of the standard concave expres¬ 
sion for electrostatic free energies can be rendered convex 
by Legendre transfornf®. Several conventions exist in the 
literature-^!. In particular we define the transform of a 
convex function c{x ) as 

g(s) = sup(sx — c(x )) (8) 

X 


We will also use the notation g(s ) = £[c](s). The trans¬ 
formation is an involution, that is the double transform 
of a convex function is the identity. 

For this paper we will be interested in the second order 
expansion of free energies about an equilibrium point. 
For this we will use an important identity linking the 
second derivatives of the functions g and c. 


d 2 c d 2 g 
dx 2 ds 2 


( 9 ) 


valid at the corresponding pair of variables (x, s ). We also 
note that the Legendre transform of a scaled function: 
ac(jx) is given by ag(s/(a'y)). 


DUAL FORMULATIONS OF THE POISSON-BOLTZMANN 
EQUATION 

In previous work^we have detailed how to pass from 
the Poisson-Boltzmann free energy expressed in terms 
of the electrostatic potential, (f> to a dual formulation 
expressed in terms of the electric displacement D. In 
particular we found 

Fdual Ffi e id T Fi ons (79) 

where 

Ffield £[eE 2 /2] = ^ (11) 


with E = —grad </>. and 

Fi orls = C[2Xk B T cosh(e/?0)](div D - p e ) (12) 
The transform of the hyperbolic cosine is easily found 
g(s) = £[2cosh](s) = ssinh _1 (s/2) — \A + s 2 (13) 
thus 

Fdual = j + k B TXg ((div D - p e ) /eX)j dr (14) 

This form as a sum of two Legendre transforms will turn 
out to be very important for a number of determinant 
identities that we derive below. We will need also the 
expansion of g(s) = —2 + s 2 /4 — s 4 /192 + .... 


as is shown in Appendix II. 
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This transformation gives a convex free energy which 
can be used in a Monte Carlo or molecular dynamics sim¬ 
ulation. It has been constructed so that the minimum of 
eq. (141 is identical to the maximum of eq. <§• The po¬ 
tential and dual formulations are linked via the equation: 


divergence operator d r i which acts on the link variable 
A . When discretizing to a simple cubic lattice there are 
3 N variables D; and N values of the divergence which we 
associate with the lattice points. We discretize the non¬ 
linear contribution to the free energy eq. (141 as follows: 


div D — p e = 2Ae sinh (/? e(f>) (15) 

which can be understood as showing that the ionic 
charge concentration in the Poisson-Boltzmann equation 
is 2Ae sinh Pecj>. 


Fions = Y g (Ez Az A - Pe) (19) 

r 

Take a second derivative with respect to the link variables 
D p and D q to find the matrix 


\q = Y d rpd rq g" (Ei AzA - Pe) (20) 

Parameterization and Scalarizing r 


It is clear that the use of a vector functional rather 
than the usual scalar form has increased the number of 
variables that need to be simulated or optimized!^ 5 . We 
show here that given a functional eq. (141 of a vector field 


we can find a related functional of a scalar field. 

When we consider the stationary point of the energy 
eq. (141 we find that 


D 

e 


k B T 


grad g = 0 


(16) 


We then identify the function k B Tg '/e = —cj> at the sta¬ 
tionary point. We see that we can make such a substitu¬ 
tion even away from the minimum to find a more restric¬ 
tive functional which is to be optimized over a smaller 
subspace: 


Fduairf = J — \-k B T\g((div egrad (j>+p e )/e A)) dr 

(17) 

Since the space is more restrictive the minimum of this 
functional is clearly greater than or equal to the func¬ 
tional expressed in terms of D, however the absolute min¬ 
imum is compatible with the parametrization, implying 
that this functional has the correct minimum free energy. 
Indeed we can go further and take linear combinations 


where g" is the second derivative of g, evaluated at r. 
Create a matrix with g" on the diagonal then eq. (201 is 
a conventional matrix product: 


= ^2 d pr(9r^rr')d r 


( 21 ) 


We identify (—9j r ) with the gradient operator since div 
and (—grad ) are mutually adjoint. 

We similarly discretize the gradient contribution to the 
free energy eq. ^ 

$r ) (C' <P'r ) ( 22 ) 

rr',1 

and take the second derivative with respect to tj> r to find 
the matrix 


Kr' = Y cz«, (23) 

l 

When we turn e/ into a diagonal matrix we find the ma¬ 
trix product: 


A rr' = (24) 

IV 


F comb = (m + 1 )F dua i^ + mF'f, (18) 

for m > 0 which all have the same, correct minimum. 

Determinant identities 


which is the discrete version of the operator 
(—div egrad ). 

We now consider the relation between the one-loop free 
energy, eq. and its naive equivalent within the dual 
formulation found by taking the second functional deriva¬ 
tive of eq. © which we write as the discretized form: 


In this section we present some matrix identities which 
link the functional determinants of operators with their 
dual equivalent. These identities are particularly inter¬ 
esting because they link operators (at least when discre- 
tised) of dimensions N x TV which are expressed in terms 
of the potential, to operators of dimensions 3 N x 3 N for 
the electric displacement. 

The identities are easiest to derive from a discretized 
form of the free energy. We will work with the discrete 


F, 


dual,loop 


k B T 

2 


Sw 

ez 


+ d T g"d 


(25) 


We will show that modulo some trivial local contributions 


to the free energy the expression eq. (251 contains the 
same physics as eq. © which we write as 


Floop — 


k B T 


ded T + d‘ 


(26) 
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where c" is the second derivative of the entropic contri¬ 
bution to the free energy. In the case of the symmetric 
electrolyte c = 2k B T\ cosh(e/3</>) 

To link these two determinants we will make use of 
the singular value decomposition theorem which states 
that a rectangular matrix, A of dimensions mxn can be 
expressed as a product 


A = UYiV* 


(27) 


where U is an unitary matrix of dimensions m x m such 
that UU* = U*U = 1; V is also unitary of dimensions 
n x n and E is diagonal rectangular with non-negative 
elements on the diagonal. 

We consider the determinant for a rectangular matrix 
A: 


\AA 7 +1 


(28) 


With the help of the singular value decomposition we can 
write 


For the specific case of a symmetric electrolyte d 2 g/ds 2 = 

l/y/A + 1?. 

As noted in a recent paper the one-loop corrected 
Poisson-Boltzmann equation contains many interesting 
physical effects including image charges in the case of 
dielectric continuitie d 21 * 26 ! Thus having a form of the 
Poisson-Boltzmann that can be simulated and that in¬ 
cludes such terms could be particularly interesting in ap¬ 
plications near surfaces. 


GENERAL SCALAR FUNCTIONALS 

Entirely different approaches to convexification have 
also been proposed in a recent series of papers 13 15 . In 
particular it was shown that it was possible to develop 
a local, convex functional which is expressed in terms of 
the scalar potential, rather than the vector field D. The 
authors give several expressions but in particular find a 
local minimizing principle with the functional 


1 + AA t \ =\1 + UT,T, t U*\ = 1 + EE T 

(29) 

r r 

= 1 + E t E = |l + a t a\ 

(30) 

Fsc= j 


W) 5 


— 2k bTX cosh (/3e<f>) 


The identity is between two matrices of different dimen¬ 
sions: to x m and n x n. We apply this theorem to 


+2A e(j> sinh(/3e0) 
, k B T 


(36) 


A = 


sinh ((div egrad </> + p e )/2e\) x 


(31) 


and find 


-Yded T ^Y + l 

— 


(32) f r 

yfd 1 y/c" 


c" 

F sc = const + / 


(div egrad <f> + p e — 2e\ sinh(/3e</>)) dr (37) 
When expanded to second order in cf> this gives: 
r3e(V^) 2 


2 / 2 


We can pull the diagonal matrices c" and e out of the 
main determinants to find 


ded 1 


d T ^d+- 
c" e 


A/3e> 

(div egrad + p e ) 2 
2e 2 /3A 


- Pe<l> 


dr (38) 


(33) 


We now make use of the relation between the curvatures 
of a function and its Legendre transform eq. |9]): g" = 
1/c", so that 


At least to second order this looks like a linear combi¬ 
nation of eq. and eq. |3|. Far from any sources in 
a background and with uniform dielectric properties we 
can re-write eq. (|38|) in Fourier space: 


ded T + d‘ 


= c 


d T g"d + - 
e 


(34) 


F = 


£( 


3 eq 2 


+ A/?e 2 4- 


eV 

2Xe 2 /3 


)l 09 


(39) 


which is our required identity. The result is rather re¬ 


markable. The dual energy eq. (141 was derived purely 


by considering stationnary properties of the free energy. 
However, if we use this in a simulation at finite temper¬ 
ature we generate automatically the correct one loop free 
energy , modulo a trivial shift in the zero. Of course we 
do not have any guarantee that any higher order con¬ 
tributions are correct - in particular in the strong cou¬ 
pling limitP^HUl the addition of terms valid in a weak 
coupling expansion may not help in correctly describing 
the physics. We can also write the dual contribution to 
the free energy as 


The question that arises at this moment is whether 
such an expression also gives rise to the correct one-loop 
potential. From the dispersion relation we deduce that 
the fluctuation potential can be written in the form 


Ffiuc = k ^~ Y log (1 


3 (<?/k) 2 + ( q/n) 4 


(40) 


which corresponds to eq. (541 in Debye-Hiickel theory. As 


„ _ k B T 

-F dual,loop — ^ lOg 


— grad g" div + - 
e 


( 35 ) 


can be expected this gives rise to a divergent integral in 
three dimensions. However it is not clear how to separate 
this expression into a Born energy with a Debye-Hiickel 
fluct uat ion potential of the form given in Appendix II, 
eq. (|58|. We must regularize by subtracting off log(g 4 ) 
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but this still leaves a Born-like energy which has a radius 
which is itself a function of the Debye length. Our con¬ 
clusion is that the scalar functionals are indeed correct 
for the ground state, but do not correspond to the correct 
excitation spectrum of the true physical system. These 
functionals can be used in a Car-Parrinello scheme where 
they are thermostated to zero temperature to avoid fluc¬ 
tuations. 


which implies that on a fine grid S ~ (1/k/i) 4 For eq. (Tt| 
the stiffness is given by S ~ L 2 /k 2 /i 4 where L is the sys¬ 
tem size, The presence of higher derivatives in a func¬ 
tional can imply a slower algorithm, even though fewer 
degrees of freedom must be integrated. 

We now move on to the second major topic of this pa¬ 
per which is the performance of optimization algorithms 
that look for minima in a free energy. 


SPECTRUM AND STIFFNESS 


SEARCHING FOR SADDLE POINTS 


In a uniform dielectric background far from charges the 
determinant eq. © can be diagonalised using a Fourier 
decomposition. The eigenvalues of the operator in eq. |5| 
are then 


0 J q = e(q 2 + k 2 ) (41) 


We now consider the spectrum in the dual picture by 
expanding the free energy eq. to quadratic order, to 
find the vector equivalent of the Debye-Hiickel theory. 


F vec = const + 


D 2 

~2e 


+ ksT 


(div D - 
4W 2 



dr (42) 


The spectrum for quadratic fluctuations decomposes into 
a longitudinal and transverse part. 


w; — ^ — w o(<7 2 + ft 2 ) (43) 

1 

u t =- 


We note that the functional forms of eq. (411 and eq. (43) 
are the same which explains, at least partly, the identity 
between the free energies at the one-loop level. 

Consider now discretizing to a grid of step h, finer than 
the Debye length. There is a gap in the spectrum eq. (43) 
and the ratio of the stiffest to the softest modes scales as 


(l/nh) 2 


(44) 


It is clear that this ratio can become large when using 
a very fine grid, in which case the system of equations 
describing the electrolyte can become stiff. When inte¬ 
grating a system with a second order algorithm in time 
(such as molecular dynamics) the number of time steps 
required to sample all modes in the system scales as \/~S. 
In a Monte Carlo simulation we can expect that the equi¬ 
libration of the electrolyte takes O(S) sweeps. Stiffness 
slows the convergence of both molecular dynamics and 
Monte Carlo algorithms. 

If we make a similar calculation for energy functionals 
such as eq. (391 we find that the equations are stiffer: 


The spectrum can be written in the form 


UJ = w 0 (l + 3 (q/n) 2 + (g/ft) 4 


(45) 


We start by noting that the minimization of a con¬ 
vex function is generally easy. Many algorithms can be 
proved to converge quickly to the correct point. As an 
example for general quadratic minima the conjugate gra¬ 
dient method for N variables converges in N iterations. 
If the object being minimized is a local free energy each 
evaluation of the energy is itself of complexity O(N). 
Thus we expect that the solution time for a model dis¬ 
cretized to N lattice points should scale as 0(N 2 ). We 
note that this may not be the optimal method of so¬ 
lution - sophisticated solvers of the Poisson-Boltzmann 
equation such as Aquasol are rather based on multigrid 
metliod ij 27 f 28 ( 

When searching for saddle points no such simple results 
can be found. Consider for instance a simple alternating 
scheme applied to two variables x and y with the energy 

/ = -x 2 /2 + y 2 /2 + axy 

For all a the saddle point remains at the origin. 

A simple algorithm is to alternate the two variables 
and optimized at each step giving the maps: 

x ^ ay y —► y (46) 

y —► —ax x —> x (47) 

A cycle of updates gives the product 



This matrix has eigenvalues 0, a 2 . We converge to the 
origin only when a 2 < 1. We learn that mixed optimiza¬ 
tion over coupled convex-concave spaces can fail when 
the coupling between variables becomes too strong in a 
way which is impossible for simple convex optimization. 

Such mixed saddle point optimization problems occur 
very often in electrostatic problems which are formulated 
in terms of the electrostatic potential coupled to some 
other, physical degree of freedom. In the field of molec¬ 
ular simulation we could consider the configuration of 
a macromolecule described by atomistic potentials. As 
a more definite example we will consider a model of a 
charged polyelectrolyte, in which a long polymer is de¬ 
scribed using an auxiliary field T, where d^ 2 describes 
the local monomer concentration 29 30 . We now consider 
a number of optimization techniques applied to such a 
system. The techniques that we try are 
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• Direct nested optimization between the polymer 
and electrostatic potential 

• Building a positive functional from squared deriva¬ 
tives 

• Legendre transform of the electrostatic free energy 

• Use of generalized scalar functionals for the elec¬ 
trostatics 

We evaluate the methods as a function of the discretiza¬ 
tion looking for methods that are rapid and converge eas¬ 
ily without extra derivative information being given by 
the user. 


A model for mixed optimization in electrostatics 


We work with a model developed to stu dy sp ontaneous 
assembly of single-stranded RNA viruse s' 31 ' :i2 [ The ge¬ 
netic material of the virus is modeled by a single poly¬ 
electrolyte chain, which must insert itself into a charged 
capsid. The system is supplemented by a symmetric elec¬ 
trolytic solution. The free energy density describing the 
interaction between the polymer, capsid and the elec¬ 
trolyte is 





(V</>) 2 + (f>(a — pe4> 2 ) 


— 2A k B T cosh (/? ecj)) 

+ k B T {y(V^) 2 + ^ 4 } 


dr . 


(48) 


The third line of this expression corresponds to the non¬ 
electrostatic interactions of the polyelectrolyte, er is the 
fixed charge of the capsid, p the charge on a monomer 
of the polyelectrolyte, a the monomer size and v the ex¬ 
cluded volume of the polyelectrolyte chain. This free en¬ 
ergy is concave with respect to </> and convex with re¬ 
spect to T. These two fields interact through the cou¬ 
pling pe^ip. In our studies we only consider systems 
with spherical symmetry, discretizing to N points, see 
Appendix I. 


Nested Optimization 


The stationary point of the functional F in equa¬ 
tion (48) corresponds to the maximum over 4> and the 


minimum over 4'. The most direct way to find such a 
point is to search independently for the two extrema: 
one optimizes iteratively the configuration while solving 
the electrostatic problem at each calculation of the itera¬ 
tive method. The outside loop of optimization calls upon 
the interior one until the saddle point is reached (see the 
algorithm 0 below). We call this method nested opti¬ 
mization. It is the kind of algorithm which is easy to 
implement if one has a pre-existing Poisson-Boltzmann 


solver which one wishes to use to study the relaxation of 
a bio-molecular system without reformulating the energy 
functionals. 


ALGORITHM 1. Nested Optimization Program 

input 

<t>*~ <t> 0 

$(-$o 
over >k 

over (j> for \k f ixe d <- \k 
minimize f ixed ) over cj> 

return (j> opt , F((j> opt ^fixed) 

<t> t— <j>opt 


4* fixed ^ (f) 

minimize F{(j)fi xf ,d, \k) over \k 
return ^ op t, fi^opt^ opt) 

*k 'it opt 

Output 'itoPT, (j>OPT, F(<j)OPT,'&OPT ) 

It is to be noted that this method differs from an it¬ 
erative scheme where both search loops are on the same 
level but included in a general loop. The schematic view 
of these two programs given in Fig. ([1} shows these differ¬ 
ences. The classical iterative method was implemented 
and tested but no convergence could be reached. We thus 
focus on the nested optimization method. 



FIG. 1. Schematic view of two iterative algorithms. On the 
left-hand side, the nested optimization method studied in this 
paper optimizes the pseudo free-energy over one field for each 
optimization step taken for the other held. On the right hand 
side, the other iterative method puts both optimization loops 
on the same level in a third encompassing loop. 


Squared gradient 

An alternative to the nested optimizations discused 
above is the use of a new functional based on the gra¬ 
dient of the free energjEl. At the saddle point 

5F\ _ /SF 

W ) *s P ,*sp ~ V ' 5 ’ 5 


= 0 (49) 
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We can build a new functional which is always positive 
and vanishes at the stationary point: 


Fderiv (^3 ^) 



The minimum of this functional yields the fields at equi¬ 
librium. The advantage of this formulation is that all 
fields can be treated equally by the minimization process 
which can be managed as a single global optimization. 


Generalized scalar functionals 


The next algorithm that we test is that based on 

- and work 


eq. (371. We denote U/) = 
with the free energy 

W,*) = / \k B T <+ ^ 


eA 




+2AfcgT | — cosh(^e0) + /3e<j> sinh ((3e<t>) 

— sinh -1 (£/2) sinh(/3e</>) + (£/2) sinh -1 (£/2) j 


dr 


The global minimum is found with a single optimization 
loop. 



FIG. 2. Convergence of free energy as a function of the num¬ 
ber of points in the discretization. Blue dots (•): Legen¬ 
dre transform. Green triangles: (A) squared gradient. Black 
squares: (□) nested loops. Red crosses ( x ): generalized scalar 
functionals. Simple minimizer without external derivative in¬ 
formation. Only the Legendre method reliably converges for 
all N. 


• Convergence of F when the number of points in¬ 
crease. 


Legendre Transform 


The Legendre transform of the electrostatic degrees of 
freedom is still possible in the presence of the extra poly¬ 
mer field. We find the locally defined free energy: 


F l (D,T) 



^(VT ) 2 + ^ 



■k B T\g((;) 


dr 


with £ = (<r — peT 2 — div D)/eA. The minimum over 
both T and D is found with a single optimization loop. 


Numerical Results 

The four functionals were programmed in Matlab. We 
use the function fminunc to search for the minimum for 
each case. This function uses a quasi-Newton algorithm^ 
and in particular the Broyden-Fletche-Goldfarb-Shanno 
update. The idea is to use a good, black-box minimizer to 
study how each formulation converges in time and accu¬ 
racy as the number of discretization points increase. The 
ideal is a formulation that converges to high accuracy 
without the need for detailed tuning of the algorithm for 
each application. We choose to stop iteration when the 
functional is converged to within an estimated fraction 
of 10 -12 . 

Three quantities are considered to study the perfor¬ 
mance of each method: 


• The derivative of the functional with respect to 
each field. We use the Iq-norm to test the validity 
of the simulation 

• Simulation time as a function of N. 


Variables are initialized to $ = 1 for the polyelec¬ 
trolyte field and </> = — 1 for the potential. The system 
size is 24 nm, where the charge of the capsid is a = OAe 
at R = 12 nm. The bulk concentration of monovalent 
ions is A = 10 mmol.L -1 and the water relative permit¬ 
tivity is e B = 80. This corresponds to a Debye length of 
3 nm. The parameters of the polyelectrolyte are set to 
a = 0.5 nm, v = 0.05 nm 3 , and p = 1. The evolution of 
the free energies as a function of N is depicted in Fig. <§; 
Fig. @ shows the derivatives on stopping. 

We work with values of N up to 400 which gives values 
of (1/nh), eq. (441, up to 50. The discretizations are 
thus rather stiff. The most robust method appears to be 
that based on the Legendre transform. Indeed, it is the 
method that converges and gives the smallest error when 
the number of points is large, Fig. (§, Fig. (§. We note 
that the squared gradient method gives poor results with 
even small N and the generalized scalar functional shows 
poor convergence for large N. 

The simulation time of each method is shown in 
Fig. © as a function of discretization; we use a recent 
intel-based desktop computer. As expected the mini¬ 
mization time scales quadratically with the number of 
variables in the discretization. We tried other, random, 
initialisations and found that in all cases the method 
based on the Legendre transform gives the most stable 


















FIG. 3. Color code as Fig. [2] Derivatives of the free energy 
for the virus model for different discretizations and methods 
on exit from the optimization loop. 

results, together with a time of calculation which is as 
good as the other tested algorithms. 



FIG. 4. Color code as Fig. [2] Total simulation time as a 
function of the number of points in the discretization of the 
free energy. All algorithms display a quadratic growth in time 
with N. The algorithm based on the double loop is much 
slower than the three other methods. Dash line slope 2. 


Stiffness and algorithms 

From our numerical experiments we see that naive op¬ 
timization of an expression with a saddle point leads to 
slow convergence. It is more surprising that some algo¬ 
rithms that can use convex optimization also give poor 
numerical results. One possible explanation for the differ¬ 
ences observed between the three convex methods is that 
the algorithms have very different stiffnesses when the 
number of discretization points increases. As noted above 


an algorithm based on the Legendre transform has a stiff¬ 
ness which increases as TV 2 , whereas both the squared 
gradient and generalized scalar functionals become much 
stiffer with scaling in TV 4 . This renders numerical codes 
much more susceptible to numerical round-off errors and 
also requires more careful stopping criteria. However as 
stated above our philosophy is that we are looking for 
free energies that are stable and easy to use without large 
amounts of algorithmic tuning from the user. 



FIG. 5. Derivatives of the free energy for the virus model 
for different discretizations and algorithms on exit from the 
optimization loop. Blue dots (•): Legendre transform, quasi- 
Newton. Cyan circles ( ): Legendre transform, trust-region. 
Red cross x : scalar functional, quasi-Newton. Magenta stars 
(*): scalar functional, trust-region. 

More sophisticated black-box algorithms are also avail¬ 
able. In particular if we are willing to calculate and 
program a routine which calculates the first derivative 
of the free energy with respect to each variable 35 we 
can find better results. In particular we use another 
algorithm implemented in Matlab the trust-region algo- 
rithrr iPI It avoids over-large steps thanks to the limit 
imposed by the definition of a trust-region, yet it main¬ 
tains strong convergence properties. The L 1-norm of the 
derivatives Fig. ([5]) and time used Fig. |6| are compared 
with the results obtained with the quasi-Newton algo¬ 
rithm. The change of algorithm allows one either to 
significantly reduce the computational time (Legendre 
transform method, blue and cyan points), or to better 
converge (generalized scalar functionals, red and magenta 
points. With the trust-region algorithm the lack of con¬ 
vergence of the free energy disappears). Thus, using the 
more sophisticated algorithm with derivative information 
renders the optimization more reliable or more efficient. 

CONCLUSION 

We have studied the quadratic fluctuations around 
equilibrium for several expressions for the free energy of 















9 



FIG. 6. Color code as for Fig. [5] Computational time for 
different methods. 


charged media. These fluctuations are closely related to 
dispersion energies and the Debye-Hiickel contribution to 
the free energy of electrolytes. The functional based on 
the field D correctly reproduces this free energy, even 
though the functional was originally proposed as a mini¬ 
mum principle. 

When used in molecular dynamics codes it is important 
to chose functionals which are not too stiff. We empha¬ 
sised that the ratio of the largest to smallest eigenvalue 
in the quadratic form is closely related to the number 
of time steps which are required to sample all modes. 
Again the functional based on the vector field D seems 
to display good scaling. 

We used a model of a virus to implement and com¬ 
pare the optimization of four free energy functionals. 
The Legendre-transformed functional shows faster con¬ 
vergence and a greater accuracy in our tests. We suspect 
that this is linked to the lower stiffness of the Legendre 
formulation of the free energy. 


Acknowledgment 

Work financed in part by the ANR grant FSCF. We 
wish to thank Derek Frydel for discussions and critical 
readings of the manuscript. 


APPENDIX: SPHERICAL SYMMETRY AND 
DISCRETISATION 

Our study of determinant identities brought out the 
importance of discretizations that conserve the adjoint 
relation between the divergence and gradient operators. 
Such identities are also important for our study of the 
minimization of the polyelectrolyte free energy. Non¬ 
compatible discretizations give rise to numerical errors 
which are different between the different formulations. 


We now show how to discretize within a spherical ge¬ 
ometry in such a way as to include important dualities 
between the discretized equations. The potentials <f> and 
'F are defined on the points of the discretisation, while 
V0, V'F and D are defined on intermediate links. We 
wish to conserve the fundamental identity 


J D • V(j) 


</div D + 


0D ■ dS 


(51) 


The discretized left hand side of this equation gives: 



(52) 


where we have used the simplest differencing scheme for 
the gradient operator and we define r n , n +i as the posi¬ 
tion of the centre of the link. Identifying both sides of 
the identity, a definition of the discretized divergence is 
obtained: 


div „ D = 


D 


n,n+\ I n ,n-(-1 


D ! r 2 

± -'n—l,n 1 n —l,n 


(53) 


The Laplacian operator is then defined as usual by V 2 = 
—div n grad . 


APPENDIX: ONE-LOOP FREE ENERGIES 

We give here the treatment of the one-loop correction 
in a homogeneous electrolyte. The free energy (com¬ 
pared with the free energy of an empty box) coming from 
quadratic fluctuations is 

( 54 ) 

9 

fVdglogU + K 2 /? 2 ) (55) 

( 27r ) Jo 

- + O{l/q 0 ) (56) 

Where we introduce an upper cut-off qo . The divergence 
in < 7 q can be balanced by a purely local self energy of the 
form 1/(2 q 2 ) corresponding to the Born energy l/87rea 
per atom, with a a real-space cut-off. The long-ranged 
Debye-Hiickel energy to the electrolyte is then 

Fdh k B Tn 3 

— = —nT (57) 

which comes from 

F D h = fast 1 + K ' 2 /? 2 ) - A/3e 2 / {eq 2 )) (58) 
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